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Chaotic dynamics of electric-field domains in periodically driven superlattices 
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Self-sustained time-dependent current oscillations under dc voltage bias have been observed 
in recent experiments on n-doped semiconductor superlattices with sequential resonant tunneling. 
The current oscillations are caused by the motion and recycling of the domain wall separating 
low- and high-electric-field regions of the superlattice, as the analysis of a discrete drift model 
shows and experimental evidence supports. Numerical simulation shows that different nonlinear 
dynamical regimes of the domain wall appear when an external microwave signal is superimposed 
on the dc bias and its driving frequency and driving amplitude vary. On the frequency - amplitude 
parameter plane, there are regions of entrainment and quasiperiodicity forming Arnol'd tongues. 
Chaos is demonstrated to appear at the boundaries of the tongues and in the regions where they 
overlap. Coexistence of up to four electric-field domains randomly nucleated in space is detected 
under ac+dc driving. 

PACS numbers: 73.20.Dx, 73.40.Gk, 47.52.+j 



I. INTRODUCTION 

Negative differential conductivity (NDC) in weakly- 
coupled narrow-miniband semiconductor superlattices 
(SL's) results in the formation of electric-field domains, 
which have i-bflen studied both experimentallju and 
theoretically.a~El The domains are stable if the doping 
or photoexcitation are large enough to form a stationary 
charge accumulation layer (the domain wall). The do- 
main wall moves from well to well as the bias increases 
and gives rise to the jumps (discontinuities) of the current 
in the stationary I-V characteristics. 

When the carrier density is not sufficiently large to 
form stable domains, but it is large enough for the uni- 
form field distribution to be unstable, periodic time- 
dependent oscillations of the current under fixed dc bias 
are possible. These oscillations have been observed in re- 
cent experiments on GaAs/AlAs SL's; they are damped 
for undoped photoexcited samplesQ and undamped for 
doped SL's without|-photoexcitationBcl According to a 
discrete drift modeljj the current oscillations are caused 
by a periodic motion of the domain wall over a few pe- 
riods of the SLjld, and this, is confirmed by photolu- 
minescence measurements. tMj Note that damped oscil- 
lations (with frequencies up to 20 GHz) hpje been ob- 
served experimentally by Le Person et alx-i in a pho- 
toexcited wide-miniband SL with strong interwell cou- 
pling. These oscillations have been interpreted in terms 
of dipole charge waves. We do not consider the miniband 
transport regime in the present paper. 

This situation is reminiscent of that found in bulk 
GaAs with NDC caused by the intervalley transfer of 
electrons, where the well-known Gunn oscillations medi- 
ated by high r field domain dynamics may appear under dc 
voltage bias.l 12 H 13 l There are several important differences 
between the Gunn oscillations and those observed in the 
SL's. A significant difference is that the space charge 
wave is a dipole in the case of the Gunn oscillations (a 



domain of high electric field with two charge layers: accu- 
mulation and depletion) and a charge monopole in the SL 
current oscillations (a wave of one charge accumulation 
layer, or the domain wall, connecting two electric-field 
domains). Another difference is that the Gunn waves are 
generated close to the injecting contact-whereas the do- 
main walls appear clearly inside the SLEJ. Notice finally 
that the different transport mechanisms determining the 
velocity of waves and the characteristic frequency of the 
oscillations give rise to different limitations in the per- 
formance of possible devices. While for the Gunn effect 
the oscillation frequency is limited by a parameter of the 
material (the intervalley scattering time), for the SL's it 
is determined by the tunneling time. Hence, the oscilla- 
tion frequency can be varied by tuning the growing pa- 
rameters (barrier widths, etc.) and/or the dc bias. The 
bias regions between different resonance peaks give rise 
to quite different frequencies, ranging between hundreds 
of KHz and several GHz over wide temperature ranges 
including room temperatureE 

In our previous papei0 we predicted the appearance 
of chaos by applying an additional small ac signal to 
the dc voltage bias in the regime where the SL exhibits 
self-sustained oscillations. In particular, when the ra- 
tio between the natural and the driving frequencies is 
fixed at the golden mean, chaotic current oscillations and 
strange attractors with a multifractal measure have been 
obtained liil 

Here we present a detailed numerical study of nonlin- 
ear current oscillations in SL's over the whole driving 
frequency - driving amplitude parameter plane. Based 
upon a relatively simple (but otherwise self-consistent) 
discrete drift model,El the shape of frequency-locked re- 
gions (Arnol'd tongues), chaotic and quasiperiodic re- 
gions is determined, thereby providing a global bifurca- 
tion picture. The largest Lyapunov exponents are cal- 
culated to identify chaotic solutions and to character- 
ize their fractal dimension. The critical line where the 
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Arnol'd tongues begin to overlap and chaos appears is 
found at very low values of the amplitude of the driv- 
ing force, unlike other periodically driven semiconductor 
systems which also display complex nonlinear interaction 
between internally generated oscillations aiicLaii external 
ac signal (e.g., GaAs Gunn diodes IB p-GeW). 

An additional point worth noticing is that the physi- 
cal mechanism of the NDC in the SL system (responsible 
for chaos) is based on a pure quantum-mechanical ef- 
fect (resonant tunneling) which is absent in the classical 
limit. Thus we may have-a case of quantum chaos with- 
out classical counterpartta which might be easier to de- 
tect ejtperimentally than that proposed by Jona-Lasinio 
et al£3. In both cases the self-consistent mean-field po- 
tential created by the charges is instrumental in inducing 
chaos, which is also the case when a quantum system is 
coupled to purely classical subsystems having a widely 
different time scale. E£l All these phenomena are different 
from the so-called quantum chaos J23 i.e., the behavior of 
quantum systems whose classical counterpart is chaotic. 

Recently chaotic dynamics in semiconductor superlat- 
tices with miniband transport has also been considered 
within a classical balance-equation approach.E!l In this 
case the chaotic dynamics is due to the negative effective 
mass of the electrons in the minibands for the regime 
of the Bloch oscillations under an external ac electric 
field (field self-consistency is also taken into account). In 
this completely different physical situation, the field was 
assumed to be spatially homogeneousEH which precludes 
consideration qithe spatial structure of chaos, in contrast 
with our workJli 

The paper is organized as follows. In Sec. II we de- 
scribe the physical model and its governing dynamical 
equations for a dark doped SL. In Sec. Ill we recall the 
main qualitative features of the self-sustained oscillations 
in a dc- voltage biased SL. The general behavior of the 
self-oscillating SL subjected to external ac and dc bi- 
ases is discussed in Sec. IV. Frequency-locked solutions 
(Arnol'd tongues) are analyzed and compared with those 
for the Gunn diode. The chaotic dynamics at the golden 
mean ratio between the natural and the driving frequen- 
cies is described in Sec. V. Different techniques to ana- 
lyze chaotic behavior (Poincare mapping, bifurcation dia- 
gram, Lyapunov exponents, phase plots, Fourier spectra, 
first return map, etc.) are discussed. Section VI draws 
the main conclusions of this work. In the Appendix we 
present the algorithm to calculate the Lyapunov expo- 
nents for our system and discuss the Lyapunov dimension 
calculated by the Kaplan- Yorke formula. 



II. THE PHYSICAL MODEL 



We consider a set of weakly interacting quantum wells 
(QW's) under high voltage biases, so that the electron 
states are localized in the wells and the tunneling process 
is sequential. The relaxation from excited levels to the 



ground state within each QW (~ 1 ps) is considered to be 
much faster than the tunneling process between adjacent 
QW's (~ 1 ns). Therefore, a single QW reaches a local 
equilibrium between two consecutive tunneling processes, 
and the state of the system can be characterized by values 
of the electric field £i(t), and the electron [density rii(t), 
with i = 1, . . . , N denoting the QW index. B 

The one-dimensional equations governing the dynam- 
ics of the doped SL are the Poisson equation averaged 
over one SL period I, Ampere's equation for the balance 
of current density, and the voltage bias conditions 



- (Si-Si-x) = - (m-No), 
e — + eniV{ti) = J, 

N 



(i) 

(2) 
(3) 



i=l 



Here e, e, and Njj are the average permittivity, the elec- 
tron charge, and the average doping density, respectively. 
The total current density J(t) is the sum of the displace- 
ment current and the electron flux due to sequential res- 
onant tunneling eriiv(£i). The effective electron velocity 
v(£) (proportional to the tunneling probability) exhibits 
maxima at the resonant fields for which the adjacent 
levels of neighboring QW's are alignedES (Fig. [I]). No- 
tice that the mechanism of sequential resonant tunneling 
is responsible for the NDC region in the velocity curve 
which in turn gives rise to the current instabilities in the 
SL. The voltage V(t) in (|h is the sum of a dc voltage V& 
and an ac microwave signal of relative amplitude A and 
driving frequency f d : V(t) = V b {l + As'm(2irf d t)}. 

The boundary condition at the first contact e(£% — 
£o) I (si) = Tlx — Vd = 8 allows for a small (8 -C Nd) neg- 
ative charge accumulation in the first well. The physical 
origin of 8 is that the n-doped SL is typically sandwiched 
between two n-doped layers with an excess of electrons, 



thereby forming a n + - 



diode.13 Then some charge 



will be transferred from the contact to the first QW cre- 
ating a small dipole field that will cancel the electron flow 
caused by the different concentration of electrons at each 
side of the first barrier. 

Note that the rate equation for the electron density for 
our model can be derived by differentiating (fil) and using 



~dt 



+ - [niv(£i) - m-\v(£i-\) 



(4) 



This is the equation of charge conservation under sequen- 
tial resonant tunneling between neighboring QW's. 

Finally, the current density J(t) in the external circuit 
under voltage bias condition can be obtained from the 
time-derivative of (^|) and Ampere's law (|^) in the form 



N 



. . dV e \ ^ 
^ =CV ~dt + N^ n]V ^ 



(5) 
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where cy = e/(lN) is the intrinsic SL capacitance per 
unit cross-sectional area. 

To study our equations, it is convenient to render them 
dimensionless, using characteristic physical quantities. 
As the unit of the electric field E = £/£\-2, we adopt 
£i_2 = A /(el), with A being the energy separation be- 
tween the first and the second electron subbands [the 
maximum of v(£)]. This yields the characteristic charge 
density no = e £]_2/(eZ) used to normalize the doping 
v = Nn / no and the electron density. The dimensionless 
velocity v(E) is obtained normalizing v(£) by its value at 
£ = £1-2, where it has a local maximum due to resonance 
in tunneling (see Fig. |l|). The others dimensionless quan- 
tities are defined as follows: the time r = t/t tun where 
ttun = l/v(£\-2) is the characteristic tunneling time; the 
dc bias V = I^/fi-jlJV; the ac bias amplitude a = AV; 
the driving frequency to = 2nfdttw 

Now substituting ([j]) and (@) into (||) we obtain a sys- 
tem of N equations for the electric-field profiles 



JY 



3=1 

— v(E{) [E l — Ei-i + v\ + awcos(wr), (6) 

with the boundary condition Eq = E\ — vb and initial 
conditions Ei(0) — V, Mi. 

As an example, we consider a GaAs/AlAs SL at T=5 
K with N=AQ, Z=13 nm, N D «1.15x 10 17 cm- 3 , A » 135 
meV, for which undamped timardependent oscillations of 
the current were first observed.u For these values one gets 
£i-2 ~10 5 V/cm, v w 0.1, and itun ~ 2.7 ns, and we take 
V=1.2 (corresponding to 14 «7.8 V). 

The system (||) is solved numerically by the fourth- 
order Runge-Kutta method using 4000 time steps per 
one oscillation period. We start with a uniform initial 
field profile and solve the equations for dc bias. After a 
short transient, the self-sustained oscillations set in and 
we switch on the ac part of the bias. Every time step we 
calculate the electric- field distribution over SL Ei(t) and 
the total current density J(t) [from Eq.(^)]. The main 
features of our numerical results are as follows. 



III. SELF-SUSTAINED OSCILLATIONS UNDER 
DC BIAS 

For pure dc case (a=0) the SL exhibits undamped 
time-periodic current oscillations when the doping den- 
sity is between some critical values v* and v** . In Fig. || 
the temporal behavior of the current starting from t=0 
is shown for different doping densities. Below v* « 0.066 
the electric-field distribution over the SL remains al- 
most uniform, so that the current exhibits no oscillations. 
Above v** w 0.175, after some transient period of the or- 
der of tt un N, stable stationary electric-field domains are 
formed. This results in an increase of the current and its 
saturation with time. 



For in-between values v* < v < v** and appropri- 
ate dc voltage bias, we find undamped oscillations with 
frequency slightly dependent on the doping. In Fig. ^ 
the frequency of the natural oscillations versus doping 
concentration is shown for different boundary parame- 
ters 6. One can estimate the natural frequency to be 
approximately equal to the inverse total tunneling time 
V(ituiiAf), which gives, after substitution of our value for 
itun ~ 2.7 ns and N = 40, a frequency about 10 MHz ip. 
close agreement with the value observed in experiment .u 
This estimation improves as N increases. All the results 
presented below were obtained for the value of 6 = 10~ 3 . 

The spatial structure of the oscillations can be seen 
from the electron density and the electric-field distribu- 
tions of Fig. ^. Starting from uniform distribution at 
time r=0, when the dc bias is turned on, after a transient 
period r « 3, a charge accumulation wave (monopole) is 
created and then it moves toward the correspondent con- 
tact. Depending on the applied voltage, it may or may 
not reach the end of the SL before it disappears and a 
new monopole is formed starting a new period of the os- 
cillation. Simulations clearly show monopole recycling 
with two monopoles coexisting during some part of one 
current oscillation period [see Fig. |(b)] . The period of 
the oscillation is mainly determined by the total travel- 
ling time of the monopole across the SL t tun N , as was 
shown in Fig. |^. 

The total current density J(t) determines the average 
state of the system (see Eq. |5|), nevertheless the infor- 
mation on the spatial dynamics is presented in the time 
dependences of the current. In the curves of Fig. || small 
current spikes can be seen during the initial stage of the 
domain formation, which are correspondent to the well- 
to- well jumping of the charge accumulation layer. The 
similar, but less clear, fine structure can be resolved in 
the periodic part of the current oscillations (become more 
pronounced after taking the derivative dJ/dt). From the 
number of spikes per period (~6 for our curves) it is 
possible to estimate the number of QW's the monopole 
moves across. 

The existence of the threshold value for the doping v* , 
above which current instabilities take place, can be un- 
derstood from the charge conservation law (^) rewritten 
in the continuum limit as 



(7) 



Taking spatial derivative and using the Poisson law for 
H one gets 



dn dn e 
dt dx e 



+ v— + -nv' E (n - N D ) = 0, 



(8) 



where the sign of the velocity derivative over the field 
v' E is crucial. In the NDC region where v' E < the last 
term in (||) gives rise to an exponential growth of the 
charge with the characteristic time r e ~ e/ (eN E \v' E \). 
For v' E > the same term gives charge relaxation to 
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quasineutrality n w Njj. Taking into account the charac- 
teristic transit time r t ~ IN /v from the convection term 
of (|sj) we can get significant growth of charge if r t > r e , 
that is NjjN > ev / (el\v' E \) . This condition is .similar 
to the critical ni-product for the Gunn effect .113 The 
treatment above is of course too rough, but it clarifies 
the influence of different parameters on the instability. 
Rewritten in our dimensionless units as 



vN > 



(9) 



this inequality implies: (i) for fixed number of the SL 
periods N there should exist the critical value for the 
doping v* above which instabilities of the current may 
occur, and (ii) for any fixed doping v one can bring the 
system to instability by increasing the length of the SL AL 
Precise bounds for v* and v** may be found elsewhere.E2l 

Notice the importance of the NDC region in the veloc- 
ity curve v(E) and the positive sign of ^ at the bound- 
ary, which are both necessary for monopole recycling. 
Besides the oscillations of the domain wall due to the 
travelling charge wave, the field in the low- and high- 
field domains also oscillates [see Fig. f|(a)]. The field 
behind the monopole is uniform in space, up to a small 
correction of the order of S near the boundary. When 
the current reaches its maximum then the field behind 
the monopole takes values on the NDC region, and those 
corrections increase exponentially in time, nucleating a 
new monopole.E-3 

The condition (Q) in the strong limit i/N 3> v/|v^| im- 
plies a large separation between the characteristic times 
Tt and r e . This can be exploited to perform an asymp- 
totic analysis of the current oscillations, which gives a 
reasonable approximation to .-the numerical simulations 
for long SL's with N > 100.EJ A key ingredient of the 
analysis is that the domain walls become shock waves 
(i.e., discontinuities moving towards the right and sepa- 
rating quasineutral regions of low and high electric field), 
whose velocity obeys an equal-area rule.E3 



IV. FREQUENCY-LOCKING UNDER DC AND 
AC BIAS 



Self-oscillating systems that are forced by an external 
oscillating signal represent an important class of coupled 
oscillators. An inherent feature of periodically forced 
nonlinear system is that the actual oscillation frequency 
depends on the amplitude of the forcing. Therefore both 
the frequency fd and the amplitude of the driving a can 
be used as control parameters to study nonlinear dynam- 
ics of our system. 

For relatively low amplitudes we could expect two pos- 
sible effects of the forcing according to the ratio between 
the driving frequency fd and the natural frequency fo- 
If this ratio is a rational number the regime of entrain- 
ment or mode-locking will be observed, if it is irrational 



the regime of quasiperiodicity-will occur because the fre- 
quencies are incommensurate.c3 

The statement above is confirmed by our numerical 
calculations. The phase diagram of Fig. ^ shows the dis- 
tribution of frequency-locked solutions (Arnol'd tongues) 
over the frequency - amplitude parameter plane. Within 
each tongue, the current is a periodic function of time 
whose actual frequency f s does not coincide with fo nor 
with fd in general. Instead it is usually related to one 
of their harmonics. Frequency-locking can be considered 
as the trend of the driven system to keep its frequency 
close to that of the unforced oscillation, fo, by taking 
the value of the harmonic of fd which is closest to fo. 
The periods of the actual oscillations are shown by num- 
bers in Fig. [| Notice that the frequency-locked solutions 
are interspersed with regions of quasiperiodicity where 
both frequencies fo and fd coexist. At a=0 the Arnol'd 
tongues arise in intervals around rational ratios of fd/ fa, 
and here f s — f . As a > increases, the tongues open 
up and eventually intersect. Then more complicated dy- 
namic behaviors (e.g., chaos) usually occur as the solu- 
tion jumps between the various overlapping modes in an 
erratic way. 

Frequency locking is a typical nonlinear phenomenon 
that can-, be found in different periodically driven 
systems. c3 In our model the Arnol'd tongues begin to 
overlap already at very low driving amplitudes (a ~ 
0.01). This is quite different of -the case of the Gunn 
diode studied by Mosekilde et aZlij In addition, we ob- 
serve 'the red shift' in the dependence of the frequency 
f s with the amplitude a, which manifests itself in the 
bending of the locked regions. The last phenomenon is 
also absent in the bulk Gunn effect and it might be a 
consequence of the delay of the natural oscillations when 
the forcing is increased. We have a qualitative argument 
to understand why the overlapping of Arnol'd tongues 
(and therefore the critical line for chaos) occurs at rel- 
atively low ac-voltages. Notice that our system has im- 
portant voltage and time scales due to its discreteness: 
recall that there are spikes superimposed to the oscilla- 
tion of the current tracking the jumps of the domain wall 
from one QW to the next one. It is plausible that over- 
lapping of Arnol'd tongues occurs when the ac voltage 
amplitude is equal to the time-averaged voltage drop per 
period, Vb/N, needed to balance the typical discrete scale 
of voltage. The latter is equal to the energy separation 
between subbands divided by the charge of the electron, 
£i-2l — A/e.c3 Unlike the ac-driven bulk semiconduc- 
tors, one could thus expect overlapping of resonances for 
a SL driven with an ac amplitude VbA ~ A/e ~ Vb/N, 
which is much smaller than the dc bias Vb- This value 
seems to give a reasonable estimation of the lowest criti- 
cal line for chaos as shown in Fig. ||. 

Let us now further describe the Arnol'd tongue dia- 
gram of Fig. |. Between the main tongues correspon- 
dent to n : 1 locking, smaller n : m intervals are found, 
representing more complex entrainments. For instance, 
between the 3 : 1 and 4 : 1 tongues one can see 10 : 3, 
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7 : 2, and 11 : 3 ratios, between the 4 : 1 and 5 : 1 
tongues there are 13 : 3, 9 : 2, and 14 : 3 ratios, and so 
on. Scanning over both parameters a and fd to analyze 
the system of 40 equations was very time-consuming, be- 
cause we calculated the largest Lyapunov exponent (see 
below) to distinguish between quasiperiodic and chaotic 
solutions. Thus we restricted ourself to scanning with a 
frequency step Afd = 0.05/o: which allowed us to detect 
the tongues n : to with at most m=3. Higher order lock- 
ing and much smaller new chaotic regions at the bound- 
aries of the tongues are expected to appear if a smaller 
step Afd is used. An additional important result in the 
mode-locking diagram of Fig. is that the richest struc- 
ture of different mode mixing and chaos is concentrated 
between the 2 : 1 and 3 : 1 tongues. In the next section 
we consider in more detail the results obtained when the 
driving frequency is fixed within that region, namely, at 
the inverse golden mean ratio. 



V. CHAOTIC DYNAMICS AT THE GOLDEN 
MEAN RATIO 



A. Bifurcation diagram 

Here we shall consider the driving amplitude a as the 
control parameter, fix the driving frequency as fd = G/o 
(where G WV5 + l)/2 S3 1.61803... is the inverse golden 
mean ratioEj), and calculate the current J as a function 
of time. 

To detect and visualize the chaotic regions in param- 
eter space, we need to define a Poincare mapping. The 
current is a good measure of the amplitude (norm) of 
the solutions, which is illustrated by the use of current 
versus voltage characteristics as bifurcation diagramsEi 
Let Td = 1/fd be the driving period. We adopt as our 
Poincare mapping (for each value of a) the current at 
times r m = mTd, m — 0,1,..., (after smiting enough 
time for the transients to have decayed) .EJ The result is 
the bifurcation diagram in Fig. 0(a) which is constructed 
as follows. For each a we compute J m = J(mTd) until 
the solution becomes periodic within a 10 -5 accuracy. 
At that time, we stop the simulation and depict all the 
J rn corresponding to one period of the solution. If the 
solution is not periodic, we eliminate the first 500 tran- 
sient points J m and depict the next 200 points. Thus, 
aperiodic solutions can be very easily distinguished from 
periodic ones in Fig. 0(a) by the large number of points 
(periods) corresponding to each value of a. Notice that 
period-2 orbits span the widest parameter region from 
a S3 0.01 up to a S3 0.085. Period-doubling cascades can 
be seen, which points out the existence of chaos near their 
accumulation points. 

The next important consideration is to discriminate 
chaotic from quasiperiodic regions, both distinguished by 
having a large number of points in the bifurcation dia- 
gram. This can be achieved by computing the largest 



Lyapunov exponent Ai. For chaotic regions Ai > 0, 
which indicates exponential divergence of nearby trajec- 
tories. For periodic solutions of our nonautonomous sys- 
tem Ai < 0, while for quasiperiodic solutions Ai = 0. 

In Fig. 0(b) we present the first and the second Lya- 
punov exponents calculated by the algorithm explained 
in the Appendix. We see that the system starts being 
quasiperiodic (Ai = 0) at a = 0, as it should be for an 
irrational frequency ratio. Then at a S3 0.005 it locks to a 
period-5 orbit (Ai < 0) terminated by some chaotic win- 
dows at a « 0.01 (Ai > 0). Notice the first appearance of 
chaos at relatively small driving amplitude (~ 1% of Vb). 
After a > 0.085 several chaotic windows can be seen, and 
then the solution becomes again quasiperiodic (as it was 
at a — 0) before locking to the driving frequency fd at 
a S3 0.145. 

More insight into the transition between chaotic and 
nonchaotic regions of the bifurcation diagram can be ob- 
tained by sweeping- up (or sweeping-down) through it, 
using the electric-field profile stored at the end of each 
simulation (for a — eio) as an initial condition for the next 
simulation (for a — ao + Aa). This approach is close to 
the process of experimental sweeping-up measurements. 
By this technique we investigated some critical points 
in the bifurcation diagram and found different pictures 
when sweep-up and sweep-down runs were made, thus 
demonstrating hysteresis. For example, such hysteretic 
behavior we obtained for the transition point between 
pcriod-2 orbit and chaos at a S3 0.085, which points out 
to the existence of a subcritical bifurcation. 



B. Current-voltage phase plots and Fourier spectra 

Fig. shows I-V phase plots and the corresponding 
Fourier spectra (FS) for some specific values of a. These 
phase plots could be measured experimentally by de- 
picting the current in the external circuit as a function 
of the instantaneous value of the ac + dc voltage bias. 
We observe different types of solutions as shown in Fig. 
7. Periodic orbits appear as simple closed loops [Fig. 
7(c),(f)]; they have few frequency peaks in the FS, cor- 
responding to the driving frequency fd and its harmon- 
ics. The quasiperiodic orbits look more complicated [Fig. 
0(a) ,(e)], but their FS are still simple. For the case of 
Fig. 0(e) (strong driving) the main peaks appear at mfd 
(to = 1,2,...). Additional small double peaks around 
(to — (with to = 1,2,...; the separation between 

the peaks of each pair is always the same for all double 
peaks) are related to one of the harmonics of the natural 
frequency /o, thereby providing coexistence of the two 
frequencies /o and fd- For the quasiperiodic case of Fig. 
0(a) (weak driving) the main peaks are at to/o with the 
same type of the double peaks, which are already related 
to fd- Thus the natural and driving frequencies exchange 
their roles depending on the strength of the forcing. Fi- 
nally, for the chaotic solutions [Fig. 0(b), (d)] the FS be- 
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come very irregular with a large number of peaks, which 
could be considered as an additional method to detect 
chaos. The FS for chaotic solutions should not be nec- 
essarily continuous. In our system the sharp frequency 
components are also present in, the FS as in the case of 
the familiar Rossler attractor .E] 

The Poincare mapping used to obtain the bifurcation 
diagram of Fig. ||(a) can be understood from Fig. [?] as the 
successive crossing of the orbit through the line V=1.2, 
where ac part of the voltage crosses zero (at upper val- 
ues of the current corresponding to the increasing volt- 
ages). For aperiodic solutions those crossing points are 
distributed over some interval (or intervals) of the cur- 
rent. 



C. First return map 

One of the main quantities observable in experiment 
is the current density in the external circuit J{t). By 
sampling the current J(t) each driving period T^, one 
obtains the data set J m . This set can be analyzed by 
means of the first return map plotting J m +i as a function 
of J m . After some transient time the resultant attractor 
for a n-period solution will be just the n separate points 
(O-dimensional object), while for an aperiodic (chaotic or 
quasiperiodic) solution it will be represented by a higher 
dimensional object. 

A closed smooth loop with regular distribution of the 
points indicates quasiperiodicity (Fig. ||a). The chaotic 
attractor is a layered (sometimes folding) structure with 
varying density of the points on different regions (Fig. 
||b), and it cap be characterized by the multifractal di- 
mension D q E3 Sometimes the chaotic attractor contains 
several separate branches almost continuously filled as 
in Fig. ||c. In this particular case the attractor has five 
branches, indicating that the solution is close to period- 
5 locking for this value of the control parameter a. The 
chaotic solution has the same symmetry as the frequency- 
locked solution which lies nearby. 



D. Spatiotemporal aspects of chaos 

So far we have only characterized the temporal aspects 
of our chaotic current oscillations leaving aside the spa- 
tial dependence of electric field and charge density. This 
dependence is a very characteristic feature of our sys- 
tem because the oscillations of the current are due to the 
dynamics of nonlinear travelling charge waves. Thus we 
analyze here what are the field and charge density profiles 
in the dynamical regimes of interest described previously. 

Under pure dc voltage bias, the SL exhibits time- 
periodic current oscillations accompanied by a periodical 
recycling of the monopolc charge wave (domain wall) in 
space as described in Sec. [II. When an ac signal is super- 



imposed on the dc voltage bias, the current can become 



chaotic for particular values of the control parameters. 
For that case the motion of the charge waves becomes 
chaotic too, as shown in Fig. ||(a). The ac part of the 
voltage causes the electric field in different parts of the 
SL to take on values in the NDC region at different in- 
stants of time. This results in irregular amplification or 
damping of the charge disturbance at different QW's. 

There are two significant new features for a SL under 
ac and dc voltage bias as compared to the pure dc case. 
The first distinctive feature is that the monopoles may 
be generated closer to the beginning of the SL, thereby 
leaving more room to their motion towards the end of 
the SL; see Fig. ^|(a). Then the peak-to- valley ratio in 
the current oscillations increases as it can be appreciated 
in Fig. 0. Amplification of the ac signal is not linked 
to chaos since it can also be observed in cases where 
there is frequency locking and an uniform electric field 
profile inside the SLE3 In fact, the amplification is most 
pronounced at the instability threshold for the doping v* , 
where the current oscillations appear. Then the electric 
field is almost uniform all the time. 

The second new feature is the qualitative change of 
the traveling wave picture for longer SL. For short SL 
(N < 80), at most two domain walls (separating three do- 
mains) can coexist at a given time (see Fig. ^(a)), whereas 
up to three coexisting domain walls (separating four do- 
mains) can be observed for long SL's during certain time 
intervals (for the case of Fig. |9|(b) at r 41.1; 46.5). In 
Fig. [lO|(a)-(c) we show the presence of one, two or three 
electric charge monopoles at different times for a 200-well 
SL under ac and dc voltage bias. We can understand this 
as follows. The width of a domain wall («6 wells for our 
particular set of parameters) is determined by the veloc- 
ity profile and it is approximately independent of the SL 
length.B Then the monopole has more room to move on a 
longer SL (under pure dc bias, monopole recycling takes 
place on about 12 periods of a 40- well SL as compared 
to about 90 periods of a 200-well SL). The disturbance 
caused by the ac bias on monopole motion is thus much 
larger for a longer SL, which results in quite different 
spatiotemporal structures under ac and dc bias [compare 
the charge densities of Figs. ||(a) and (b)]. For a short 
SL the chaotic behavior can be associated with a chaotic 
domain wall dynamics that resembles the dc voltage bias 
case: most of the time there is only one domain wall 
which moves to the end of the SL and disappears, and 
about that time another monopole is generated, some- 
times quite close to the beginning of the SL [see Fig. 
0(a)]. On the other hand, for a long SL the graph of 
the electron density is disjoint: In addition to long-living 
waves traveling over almost the whole SL, there are short- 
living waves existing only at the beginning of the SL. The 
two types of waves are distributed chaotically in space. 

We have not found more than three monopoles by in- 
creasing further the SL length N. This may be due to 
the fact that there is not enough charge in the SL (de- 
termined by doping) to provide more than three jumps 
in the electric field. It might be possible to observe more 
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complicated field structures in parameter regions where 
the current instabilities exist at much higher electron 
densities. 

The spatially chaotic nature of the solutions can be 
illustrated by picking two far-away QW's and depicting 
the simultaneous values of the electric field at them after 
each period of the driving force Tj. The resultant attrac- 
tors will be very similar to those obtain ed fo r the first 
returnjiiap J m +i — J m described in Sec. VC (cf. Fig. 2 
in Ref Jla with Fig. || of the present paper) . 



VI. CONCLUSIONS 

Dynamic properties of the high-field transport in 
weakly-coupled superlattices under ac and dc biases 
has been studied numerically within the simple self- 
consistent discrete drift model. A nonlinear interaction 
between an internally generated periodical motion of the 
accumulated charge wave and an external microwave sig- 
nal gives rise to a frequency-locking, quasiperiodicity or 
chaos depending on the external driving parameters. The 
calculated phase diagram of the frequency-locked solu- 
tions (Arnol'd tongues) shows the following new features: 
(i) the first overlapping of tongues giving transition to 
chaos occurs under very weak driving (a ~ 0.01); (ii) 
there is a bending of tongues under strong driving ('red 
shift' in frequency) which we associate with the delay in 
natural oscillation frequency, and that is a consequence 
of our spatially extended system. 

Besides the chaotic dynamics, the microwave forcing 
was found to give an amplification of the self-sustained 
current oscillations, which is most pronounced at the in- 
stability threshold for the doping v* , where the oscilla- 
tions appear. 

It would be of interest to verify our predictions by 
making measurements in currently available n-doped 
GaAs/AlAs samples forming n + -n-n + diodes. These 
samples exhibit self-sustained oscillations under pure dc 
voltage bias,B and are thus suitable candidates for ob- 
serving chaos when an appropriate microwave signal 
is superimposed on the dc voltage bias. More highly 
doped samples present multistable stationary solution 
branches, (corresponding to coexisting static electric field 
domainsa) in their current - voltage characteristics under 
dc voltage bias. Studying the response of these samples 
to ac and dc voltage bias would also be of interest. Pre- 
liminary calculations show rich spatiotemporal behavior 
including chaos. 
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APPENDIX: LYAPUNOV EXPONENTS 

Consider the ^-dimensional phase space for the 
electric-field vector Ei. The evolution of its trajectory 
Ei(r) is described by the set of nonlinear equations (||). 
Since we want to know how trajectories infinitesimally 
close to the real trajectory of our system diverge from it, 
we can use a linearized version of the equations m) : 



dej 
dr 



N 



I J2W(Ej) „J e 3 + v(£J)(e, - e^)] 
i=i 

v'(El)nl ei-viETXei-e^t). 



N 



(Al) 



Here ii, with i = 1, . . . , N, is the vector of the distur- 
bances of the electric field. The velocity v(E[), its deriva- 
tive v'(EJ), and the electron density nj = EJ — EJ_ 1 + v 
are calculated on the fiducial trajectory Ei(r). The 
boundary condition for the system (Al) is eo = e\. 

To calculate the largest Lyapunov exponent Ai is very 
easy. We take an arbitrary initial vector ej(0). Then 
both the nonlinear system (^|) and the linearized system 
(Al) are advanced forward in time simultaneously. 

Denoting by ||ej(r)|| the euclidean norm of the vector 
e, at time r, the largest Lyapunov exponent is given by 



, 1, \Mt)\\ 

Ai = lim — in T7— — t- 77. 

r^oor ||e;(0)|| 



(A2) 



The linear system without nonlinearities will give an 
exponential growth (relaxation) of the solution for a pos- 
itive (negative) Ai. Therefore, the perturbations will 
need to be renormalized from time to time to prevent 
overflow (underflow), because they are only represented 
with a finite floating point numbers in the computer. To 
avoid this we use the algorithm of Benettin et al£3 and 
renormalize our perturbation vector after each half pe- 
riod Td/2 of the ac voltage bias. Then Ai will be given 
by 



Ai 



lim 



^ In IK- 1 

3=1 



(A3) 



where dj is the vector perturbation growth during the 
jth renormalization period. Fig. |ll| shows the results 
for Ai(ife) obtained after computing of k renormalization 
intervals. In the limit k — > oo X\ converges to its limiting 
value (positive, zero, or negative) according to the control 
parameter a. Notice a little difference between a=0.01 
and a=0.0102 gives a drastic difference in values of the 
Lyapunov exponents from negative to positive. 
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To look at our chaotic solutions from the point of view 
of their dimensionality we have to compute the next Lya- 
punov numbers. Of course, we don't need all the Lya- 
punov spectrum of 40 numbers to characterize our sys- 
tem. Kaplan and Yorke defined a quantity—called the 
Lyapunov dimension Dl, given by a formulaEf] 



K 



D, 



K 



K+l 



(A4) 



»=i 



where K is the largest integer such that Y^f=i — ( au 
Xi are arranged in decreasing sequence). Since different 
Lyapunov exponents characterize the stretching and con- 
tracting of phase space in different directions, the Lya- 
punov dimension is then the number of vectors in phase 
space needed to describe an infinitesimal volume that 
remains constant on the average. Dr, is-related to the in- 
formation dimension of the system Dil3 and can be used 
to characterize the fractal dimension of the associated 
chaotic attractor. 

To calculate the first n exponents Ai, . . . , A„ we de- 
fine n initial perturbation vectors e^(0), which are or- 
thonormal, so that they form n-dimensional sphere. As 
the system is integrated forward in time, the sphere of 
trajectories will become an ellipsoid since different direc- 
tions in phase space will expand or contract at differ- 
ent rates. The Lyapunov exponents A; are determined 
by the rate of expansion or contraction of the principal 
axes of the ellipsoid averaged over the entire attractor. 
One additional problem arises here. Since all the vectors 
e$ (t) tend to line up in the direction of the greatest 
growth we have to use the Gram-Schmidt orthonormal- 
ization procedure to separate the vectors into orthogonal 
components. The Lyapunov exponents Aj are then cal- 



culated by the same formula ( |A3|) as for Ai, with d' 



(n)i 



being the orthonormalized perturbation growth of ||e 
during jth orthonormalization period. 

Fig. [l2| shows the temporal convergence of the first 
three Lyapunov exponents for the particular chaotic so- 
lution (a=0.09). The estimated from those values Lya- 
punov dimension Dl w 2.208 indicates low-dimensional 
chaos for chaotic dynamics in our SL. By further increas- 
ing the length of the SL N the characteristic Lyapunov 
dimension does not grow. Since we did not see a tran- 
sition from low- to high-dimensional chaos, it could be 
concluded that our system is not extensively chaotic. We 
have not observed hyperchaos either, where the second 
Lyapunov exponent is positive. 



On leave from Institute of Semiconductor Physics, National 
Academy of Sciences, Pr Nauki 45, 252650 Kiev, Ukraine. 
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FIG. 1. Dimensionless effective velocity (proportional to 
the tunneling probability) as a function of the electric field. 
The point indicates the electric-field value corresponding to 
the dc bias V = 1.2 used in the calculations. 

FIG. 2. Temporal evolution of the current J (under dc 
bias) for different values of the doping concentration v in- 
dicated at the right margin with increasing order from the 
bottom to the top. The boundary condition <5=10 -3 . 

FIG. 3. Frequency of the natural oscillations /o under dc 
bias (scaled on the total travelling time over the whole SL 
ttuuN) vs doping concentration v for TV = 40 and for different 
boundary parameters 5. Out of the regions marked by dashed 
vertical lines oscillations disappear. 
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FIG. 4. Spatiotemporal distributions of electric field (a) 
and electron density (b) for self- sustained oscillations under 
dc bias. Quantum well index is denoted by i. 

FIG. 5. Phase diagram showing the distribution of fre- 
quency-locked solutions over the driving frequency - driving 
amplitude parameter plane. Driving frequency fd is in units of 
the natural oscillation frequency fo ps 0.0235(l/ttun). Grey 
color corresponds to quasiperiodicity, lighter color with dif- 
ferent shades corresponds to periodical solutions (number of 
periods is shown by numbers), chaos is marked by black color. 
G is the inverse golden mean. 

FIG. 6. Bifurcation diagram of the current obtained by 
means of Poincare mapping each driving-force period Td (a) 
and the first two Lyapunov exponents (b), both versus the 
driving-force amplitude a for the golden-mean ratio between 
natural and driving frequencies. Windows of chaotic solutions 
are marked by arrows. Lyapunov exponents are scaled on the 
driving period Td- 

FIG. 7. Phase current-voltage plots (left) and the corre- 
sponding Fourier spectra of the current (right) for different 
driving amplitudes a: 0.005(a); 0.0102 (b); 0.05 (c); 0.09 (d); 
0.144 (e); 0.16 (f). 

FIG. 8. First return map for current density J m at differ- 
ent driving amplitudes a: 0.14(a); 0.09(b); 0.101(c), indicat- 
ing quasiperiodicity (a) and chaos (b),(c). 



FIG. 9. Chaotic propagation of the charge accumulation 
waves in a 40-period SL for a=0.09 (a) and in a 200-period 
SL for a=0.12 (b). 



FIG. 10. Nonuniform distribution of electric-field domains 
for a 200-period SL: (a) a=0 at the time moment r = 3.20 
(one monopole); (b) a=0 at r = 4.35 (two monopoles); (c) 
a=0.112 at t = 18.05 (three monopoles); 



FIG. 11. Temporal convergence of the largest Lyapunov 
exponent Ai for different driving amplitudes a: 0.01 (peri- 
odic); 0.0102 (chaotic); 0.05 (periodic); 0.1395 (quasiperi- 
odic). Time tu = kTd/2 is in units of the renormalization 
intervals. Lyapunov exponents are scaled on the driving pe- 
riod T d . 

FIG. 12. Temporal convergence of the first three Lya- 
punov exponents for chaotic solution at a=0.09. Lyapunov 
exponents are scaled on the driving period Td. 







